
#include "polyfit.h"

#include <gsl/gsl_errno.h>
#include <gsl/gsl_machine.h>
#include <gsl/gsl_math.h>
#include <math.h>

int cubic_minimizer (double a, double fa, double dfa,
                     double b, double fb, double dfb,
                     int extrapolate, double *alpha_min)
{
  double z, w;
  double denom;
  double D;
  
  if(gsl_isnan(a) || gsl_isinf(a) || gsl_isnan(fa) || gsl_isinf(fa) || gsl_isnan(dfa) || gsl_isinf(dfa) || 
     gsl_isnan(b) || gsl_isinf(b) || gsl_isnan(fb) || gsl_isinf(fb) || gsl_isnan(dfb) || gsl_isinf(dfb))
    return GSL_EDOM;
  
  z = 3.0 * (fa - fb) / (b - a) + dfa + dfb;
  D = z * z - dfa * dfb;
  
  if(extrapolate == 0)
    {
      if(D < 0.0)
	return GSL_EDOM;
    
      if(a < b)
	w = sqrt(D);
      else
	w = -sqrt(D);
    }
  else
    {
      if(a < b)
	w = sqrt(GSL_MAX(D, 0.0));
      else
	w = -sqrt(GSL_MAX(D, 0.0));
    }
  
  denom = dfb - dfa + 2.0 * w;
  
  if(fabs(denom) < GSL_DBL_EPSILON)
    return GSL_EZERODIV;
  
  *alpha_min = b - (b - a) * (dfb + w - z) / denom;
  
  return GSL_SUCCESS;
}
 
int quad_minimizer1 (double a, double fa, double dfa,
                     double b, double fb,
                     double *alpha_min)
{
  if(gsl_isnan(a) || gsl_isinf(a) || gsl_isnan(fa) || gsl_isinf(fa) || gsl_isnan(dfa) || gsl_isinf(dfa) || 
     gsl_isnan(b) || gsl_isinf(b) || gsl_isnan(fb) || gsl_isinf(fb))
    return GSL_EDOM;
  
  double denom = 2.0 * (fa - fb + (b - a) * dfa);
  
  if(fabs(denom) < GSL_DBL_EPSILON)
    return GSL_EZERODIV;
  
  *alpha_min = a + (b - a) * (b - a) * dfa / denom;
  
  return GSL_SUCCESS;
}

int quad_minimizer2 (double a, double dfa,
                     double b, double dfb,
                     double *alpha_min)
{
  if(gsl_isnan(a) || gsl_isinf(a) || gsl_isnan(dfa) || gsl_isinf(dfa) || 
     gsl_isnan(b) || gsl_isinf(b) || gsl_isnan(dfb) || gsl_isinf(dfb))
    return GSL_EDOM;
  
  double denom;
  
  denom = dfa - dfb;
  
  if(fabs(denom) < GSL_DBL_EPSILON)
    return GSL_EZERODIV;
  
  *alpha_min = b + (b - a) * dfb / denom;
  
  return GSL_SUCCESS;
}
